Classical density-functional theory of inhomogeneous water including explicit 
molecular structure and nonlinear dielectric response 
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We present an accurate free-energy functional for liquid water written in terms of a set of effective 
potential fields in which fictitious noninteracting water molecules move. The functional contains an 
exact expression of the entropy of noninteracting molecules and thus provides an ideal starting 
point for the inclusion of complex inter-molecular interactions which depend on the orientation 
of the interacting molecules. We show how an excess free-energy functional can be constructed 
to reproduce the following properties of water: the dielectric response; the experimental site-site 
correlation functions; the surface tension; the bulk modulus of the liquid and the variation of this 
modulus with pressure; the density of the liquid and the vapor phase; and liquid-vapor coexistence. 
As a demonstration, we present results for the application of this theory to the behavior of liquid 
water in a parallel plate capacitor. In particular, we make predictions for the dielectric response of 
water in the nonlinear regime, finding excellent agreement with known data. 

PACS numbers: 71.15.Mb, 05.20.Jj, 77.84.Nh, 87.15. kr 



o 

■4-* 



i 

C 

o 
o 



> 
o 

oo 

(N 
O 
On 
O 



X 



INTRODUCTION 

Water is the most important liquid on earth. Although 
its uniform phase already exhibits many fascinating prop- 
erties and is currently an area of intensive research, the 
importance of the inhomogeneous phase in the chemi- 
cal and biological sciences is even greater: as a solvent, 
water is ubiquitous in physical chemistry, biochemistry 
and electrochemistry. The associated phenomena include 
hydrophobic interactions 0, 01, protein folding 0, 0], 
the electrochemical interface [5j, and phase transitions 
of confined liquids @, 0|. Due to the complex inter- 
play of hydrogen bonding, long-range polar interactions, 
and short-range excluded volume effects, developing a 
tractable theory to describe the solvent in the aforemen- 
tioned systems remains a challenge Q. 

Despite the importance, inherent interest, and exten- 
sive experimental study of water, most existing theories 
of the inhomogeneous phase of water either lack the ac- 
curacy to describe the above phenomena in a quanti- 
tatively satisfying way or become computationally pro- 
hibitive. The latter is particularly true for applications 
to complex systems, such as dissolved biomolecules or 
electrochemical interfaces. Simple dielectric continuum 
theories 0, [l(| capture the dielectric response of the 
solvent at long length scales, but fail to describe more 
subtle effects at interfaces, such as the reorganization 
of the hydrogen-bond network or the binding of water 
molecules to solutes. Atomistic theories, on the other 
hand, describe such local effects quite well, but require 
the explicit treatment of many solvent molecules and also 
demand long simulation times to calculate statistically 
reliable thermodynamic averages. Since the full-fledged 
ab initio molecular dynamics description of the solvent is 
only feasible for the smallest solutes, generally additional 
approximations become necessary in practice: Bernholc 



and coworkers, for example, have introduced a clever 
scheme which employs the full ab initio description only 
for a few solvent molecules adjacent to the solute and 
use rigid geometries and frozen electron densities for the 
remaining molecules In other studies, combination 



of ab initio methods for the solute with classical force 
field approaches [12, lU 14 1 for the solvent yields addi- 
tional simplifications. Despite the simplifications which 
such hybrid approaches allow, the need to compute many 
configurations to sample phase-space properly ultimately 
limits the size of systems which can be studied and thus 
the applicability of such approaches. 

Ideally, one would have a rigorous description of the 
solvent which is able to combine the advantages of con- 
tinuum theories (large system sizes and implicit thermo- 
dynamic averaging) with an explicit geometric descrip- 
tion of the solvent molecule. Here, we show that just 
such a theory can be developed for liquid water by start- 
ing with the quantum mechanical free-energy functional 
for both the electrons and nuclei which comprise the liq- 
uid and, by "integrating out" the electrons, proceeding 
to construct a density functional in terms of atomic site 
densities alone. Such "classical" density-functional the- 
ories, which are founded on a number of exact theorems 
[lil ! have been applied successfully to the study of 
simple liquids in the past 17, HI]. Historically, applica- 
tion of this method begins with a hard sphere reference 
system that is usually augmented by terms that capture 
weak long-range attractive forces [19j. Unfortunately, a 
hard sphere reference system is a poor starting point for 
the description of water because of the strong anisotropic 
short-range interactions arising from the molecular struc- 
ture, including effects such as hydrogen bonding. 



To remedy this, Chandler and coworkers [2fJ, [2l|, [22 1 
introduced a density-functional theory for general molec- 
ular liquids in terms of a set of densities, one for each 
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"interaction site" on the molecule (typically atomic cen- 
ters). Although their theory successfully predicted the 
correct hydrogen-bonded structure of ice [22| , it suffered 
from a significant flaw: to make calculations tractable, 
the "united atom approximation" , which takes all of the 
sites on the molecule to sit at a single point, had to be 
invoked. This severe oversimplification of the molecular 
structure prohibited the description of dielectric proper- 
ties and stalled the further development of this promising 
approach. In 23J, we eliminated the need for the united 
atom approximation by writing the free energy of a gen- 
eral molecular liquid as a functional of a set of effective 
potentials in which fictitious non-interacting molecules 
move. This advance allows for the exact evaluation of 
the entropy associated with the geometric structure of 
the molecules in a numerically feasible way and forms 
the basis for the construction of a new class of accurate 
density functionals for fully interacting molecular liquids. 

In this work, we apply our general theory of molec- 
ular liquids to water and demonstrate how an accurate 
free-energy functional can be developed to reproduce the 
equation of state, the experimental site-site correlation 
functions, and the macroscopic dielectric response. The 
organization of the work is as follows. Section 2 gives 
a brief review of the general theory of molecular liquids 
introduced in [23|. Section 3 discusses the application 
of this approach to the particular case of water and in- 
troduces new computational techniques to simplify cal- 
culations of triatomic liquids. Section 4 presents results 
obtained using this new microscopically-informed contin- 
uum theory to study the behavior of water in a parallel 
plate capacitor. Here, we explore the differences between 
capacitor plates consisting either of purely repulsive hard 
walls or of walls with a local region of attraction. We also 
explore the dielectric response of water in both the lin- 
ear and the nonlinear regimes of the applied field. Com- 
paring our results to previous works based on explicit 
molecular simulations, we find quite encouraging agree- 
ment. Section 5 then concludes this work. 



we wrote 

M . 
a=l •* 

+U[n].(l) 

Here, </> a (r) is the site-dependent external potential (with 
a referring to the different atomic sites on each molecule) , 
fi a is a site specific chemical potential and t/[n] denotes 
the excess free energy due to inter-molecular interac- 
tions, which is a functional of the set of site densities 
n = {m(r), ...,nAf(y)}- Finally, f2( m )[>t] is the grand 
free energy of an assembly of noninteracting molecules in 
a set of effective potentials This free energy is known 
explicitly and is exactly 

n< n< >[¥] = -k B Tm J d 3M r s({r a }) e -^" i*»( r «', (2) 

where ni is the reference density at vanishing chemical 
potentials (which we will later take to be the bulk liquid 
density) and s({r Q }), which describes the geometry of 
the molecule, is the intra-molecular distribution function. 
Note that we take the densities in ([1]) to be themselves 
functionals of the effective fields via 



n Q (r) 



(3) 



DENSITY-FUNCTIONAL THEORY OF 
MOLECULAR LIQUIDS 



Finally, the effective potentials which minimize de- 
termine the equilibrium site densities and allow for the 
calculation of various equilibrium properties of the liquid. 

The introduction of effective potentials instead of site 
density fields as fundamental variables in the density- 
functional theory of classical molecular liquids allows for 
the first time a computationally tractable approach to 
treat the non-interacting system exactly. This parallels 
the introduction of Kohn-Sham orbitals as fundamental 
variables in the density-functional theory of electronic 
structure In both cases, the construction of highly 
accurate density functionals is made possible by map- 
ping a system of interacting particles onto a fictitious 
system of noninteracting particles with the same equilib- 
rium densities. We mention that our approach of mini- 
mizing the free energy with respect to a set of effective 
potentials is known in the electronic structure context as 
the optimized effective potential method pBj . 



Kohn-Sham Approach and Entropy Functional 

In [23|], we showed that the grand free energy £1 of 
an assembly of interacting identical molecules (consist- 
ing of M atoms or "interaction sites" each) can be ex- 
pressed as a functional of a set of effective potentials 
\I> = {Wi(r), $if(r)}, in which a corresponding set 
of fictitious non-interacting molecules move. Specifically, 



Excess Free Energy Functional 

In [23| , we also present a recipe for the construction of 
the excess free energy functional in (|T|). As usual with 
density-functional theories, this construction is difficult. 
As in our previous work, we follow the lead of Kohn and 
Sham and construct a free-energy functional that repro- 
duces established results for the homogeneous phase in 
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the limit of vanishing external fields. Namely, we expand 
U in a power series about the uniform liquid, and for- 
mally group all terms except the quadratic part, which 
is analogous to the Hartree energy in electronic-structure 
theory, into an unknown functional F ex [n] , which is then 
the analogue to the electronic exchange-correlation func- 
tional, whose form, being unknown, must ultimately be 
treated in some approximate way. 

To proceed, we separate the kernel of the quadratic 
part of U into a long- wavelength part K ai , known an- 
alytically for the case of rigid molecules, and a remain- 
der Cq, 7 , to be determined from experimental data. The 
next paragraph describes K ai in detail, and the set of 
functions C Q7 will be constructed below such that the 
resulting functional reproduces the experimental corre- 
lation functions in the uniform phase. When these two 
parts are determined, they combine with F ex to give U 
as 

M 

U[n] = J d 3 rj d 3 r' n a (v) {K aj (r, r') 



a, 7—1 



+C Q7 (r,r')Wr')+ J F ie >]. (4) 



The derivation of K aj [23j considers a collection of rigid 
neutral molecules with each interaction site carrying a 
partial charge q a with no assumptions regarding the form 
of the non-electrostatic part of the interaction — the the- 
ory neither is perturbative in the interaction strength 
nor presupposes a decomposition into pairwise interac- 
tion potentials. In [23| , we then prove rigorously that, 
for a certain class of liquids including all diatomic liq- 
uids and also water, the leading order term in a long- 
wavelength expansion of d^U, the Hessian of U for rigid 
molecules, is 



1 



c(ni) _ I 



— q a q^, (5) 



F ex then captures all of the internal energy of the uni- 
form phase (and the long-wavelength limit of the non- 
singular part of the Hessian d 2 U) and can be expressed 
as F ex = Vf ex (n), with V being the volume and n the 
average molecular density. Due to the presence of mul- 
tiple density fields, generalizing the above expression for 
F ex to the inhomogeneous case is more difficult than for 
the analogous exchange-correlation energy in electronic 
structure theory. Moreover, because of the strong corre- 
lations induced by excluded volume effects, purely local 
excess functionals fail to describe the liquid state [2(| . 
We therefore approximate F ex with a simplified ansatz 
in the spirit of weighted density-functional theory [26j |. 
but generalized to multiple species, 



F ex \n] = 



•E 



^& 7 n 7 (r) 



(6) 



where we introduce the weighted densities n 7 (r) = 
J dV^r 2 ,) -3 / 2 exp(— |r — r'| 2 /r 2 )n 7 (r'), with r being 
a "smoothing" parameter ultimately fit to the experi- 
mental surface tension. To reduce to the correct form in 
the uniform case, pi and & 7 must fulfill ^2 i p% = 1 and 

Y^' .>>'■ I- 

For the scalar function f ex (n) , we use a polynomial fit 
to various bulk thermodynamic properties of the liquid. 
The condition that C vanishes in the long-wavelength 
limit subsequently fixes pt and b^. For a given ro, 
this then completely specifies our approximation to F ex . 
Next, relating K + C + d 2 F ex to the density-density cor- 
relation functions through the Ornstcin-Zcrnike relation 
gives the matrix function C for a given r . Finally, we 
can determine the smoothing parameter ro by adjust- 
ment until calculations of the liquid-vapor interface give 
the correct surface tension. 



APPLICATION TO WATER 



where e is the macroscopic dielectric constant and e^ ni ' is 
the dielectric constant of a system with intra-molecular 
correlations only. Because K ai now captures the singular 
long- wavelength features of the response function, C ail 
which we define as the difference between d^U (given 
by the experimental correlation functions) and the sum 
of K ai and d^F ex , must now be smooth near the ori- 
gin and thus amenable to numerical approximations. As 
discussed below, we define F ex such that its Hessian 
matches the constant term in the long-wavelength ex- 
pansion of <9 2 [/ . Therefore, the functions C Q7 have to 
vanish for small k. 

Next, to approximate F ex , we begin by noting that in 
the case of zero external fields all densities are equal and 
the first quadratic term (K) in fH vanishes because of 
charge neutrality. Anticipating that, by definition, the 
matrix function C vanishes in the long- wavelength limit, 



Intra-molecular Free Energy 

Assuming rigid bonds of fixed length (and thus rigid 
angles as well) the intra-molecular distribution function 
of a water molecule (normalized to integrate to the vol- 
ume of the system) is 



s(r ,ri,r 2 ) 



5(\v w \-B)8(\v 20 \-B) 
AttB 2 AttB 2 
x2<5(?i • r 2 o ^ cos6 B ), 



(7) 



where B — 1 A is the oxygen-hydrogen bond length taken 
from the SPC/E [27| model of water, Ob is the bond angle 
between the hydrogens and taken to be the tetrahcdral 
angle like in the SPC/E model. Also, the label refers 
to the oxygen atom and the labels 1 and 2 to the two 



hydrogens, and we define 



Yj. Contrary to 
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the case of diatomic molecules, where the assumption 
of rigid bonds turns f2^ n ^ into a convolution, numerical 
evaluation of @ proves more challenging in the triatomic 
case because of the complex, nine-dimensional form of s. 

To address this, we observe here that fl^ m "> for water 
can also be brought into a "convolution form" by mak- 
ing use of the following mathematical identity, which is 
readily derived using the addition theorem for spherical 
harmonics, 



<5(v 



cosC) = J2 2vrP ; (cosC)^„(v)y ;m (w), (8) 



i m 



where v and w are unit vectors, £ is an arbitrary angle, 
Pi is a Legendre polynomial and the Y; m 's denote the 
spherical harmonics. (To prove ([5]), first expand the delta 
function as a sum of pairs of Legendre polynomials on 
the interval [—1; 1] and then use the spherical- harmonic 
"addition theorem" to express Pi (v • w) as a sum of pairs 
of spherical harmonics.) 

Substituting ([5]) and ([7]) into results in a computa- 
tionally efficient form for the non-interacting free energy 



Pi 


Pv 


B, 


dBt/dP 


[kg/m 3 ] 


[kg/m 3 ] 


[GPa] 


(dimensionless) 


997.1 


0.023 


2.187 


5.8 



TABLE I: Experimental inputs to construction of f ex from 
[28l ] and : liquid density (pi), vapor density (p v ), liquid 
bulk modulus (Bi) and derivative of modulus with respect to 
pressure (dBi/dP). 



fa 


[hartree x bohr -3 ] 


3.7608 5540 7008 46 x 10" 18 


h 


[hartree] 


-9.9058 0951 0335 88 x 10~ n 


h 


[hartree x bohr 3 ] 


8.6971 6797 7040 56 x 10" 4 


h 


[hartree x bohr 6 ] 


-2.5454 6838 7756 12 x 10 3 


U 


[hartree x bohr 9 ] 


8.9610 5939 0318 49 x 10 5 


h 


[hartree x bohr 12 ] 


-1.2391 1957 0697 09 x 10 s 


h 


[hartree x bohr 15 ] 


6.4173 3045 2123 21 x 10 9 



TABLE II: Coefficients of polynomial parametrization of f ex 
given to all digits used in our software in order to ensure 
numerical reproducibility of our results. Any individual coef- 
ficient contains at most two or three significant figures. 



n im) = -k B Tni J tfV-^ ^ 
x^4^(cos^)/^(r )/f (r ), (9) 



I in 



where fj^ and , given by 



AttB 2 



Yi m (r 20 )e-^^\l0) 



are now convolutions and can be efficiently evaluated 
with fast-Fouricr-transform techniques. Transformation 
to Fourier space greatly simplifies the above convolutions, 



/«(k) = (-i)'ii(|k|B)Y4(fc).F[ 



H)^(|k|-B)Y m (k)JF[e 



,-/3*i(r)l 
-/3* 2 (r)i 



(11) 



where ji denotes the spherical Bessel functions of the first 
kind and J r [exp(— /34 r i(r))] denotes the Fourier transform 
of exp(-/3*i(r)). 

To evaluate (|9]) numerically, we first choose a maximum 
value of I, Imaxi after which we truncate the infinite sum. 
In general, the choice of l max depends on the density pro- 
file that has to be represented. We find in the capacitor 
calculation presented below that the interaction site den- 
sities for water adjacent to the capacitor wall are fairly 
smooth in the linear response regime and can be suffi- 
ciently resolved with l m ax = 10. If strong external fields 
are applied, sharp features in the density profile develop 
and we use l m ax = 40 to be absolutely sure of a highly 
converged description. 



Intermolecular Free Energy 

First, we approximate f ex (n) as a sixth-order polyno- 
mial, 



(12) 



and adjust its coefficients to reproduce the seven condi- 
tions represented by (a) the thermodynamic stability of 
liquid and vapor phases, (b) their coexistence, (c) the 
experimental liquid and vapor densities, (d) the experi- 
mental bulk modulus of the liquid (E>i ) and (e) the deriva- 
tive of the bulk modulus with respect to the pressure P 
at standard temperature and pressure (T = 25 °C and 
P = 101.325 kPa = 1 atm). Table U summarizes the 
experimental input used to determine the coefficients / p , 
and Table HT1 gives the actual numerical values of the co- 
efficients used in our calculations below. (Clearly, not 
all of the digits given in the table are significant, we 
give them only to make our results numerically repro- 
ducible.) The values in the table are given in atomic 
units (1 hartree w 27.21 eV, 1 bohr w 0.5291 A) as this 
is the system which our software employs. 

Next, we can determine the coefficients in the inte- 
grand in from knowledge of the constant term in the 
long- wavelength expansion of the Hessian d^JJ . The en- 
tries in this term (which is a matrix) are related to various 
material response properties such as as the bulk modu- 
lus, but unfortunately, we do not have access to data for 
all of the material properties. In the absence of data, we 
take this matrix to be proportional to its value for the 
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non-interacting case with a proportionality constant set 
to ensure the correct bulk modulus. This then fixes the 
integrand in ^ and allows it to be written entirely in 
terms of rational numbers in the form, 



iq 3 Q 

27^ / rip + n 1 + n 2 
20 J \ 3 



(13) 



Next, to construct the matrix K we use the partial 
charges of the SPC/E water model 27| . i.e. q% = q% = 
—qo/2 — 0.4238 e (e being the proton charge) and the 
fact that the dielectric constant of the liquid with intra- 
molecular correlations only is e'™^ = (1 — inpmp 2 /3) _1 , 
where p is the dipole moment of an SPC/E water 
molecule. Also, we employ the experimental value [28j ]. 
e = 78.4, for the macroscopic dielectric function of liquid 
water at standard conditions described above. 

The use of K as leading order term is only justified 
for long-wavelength (small k) properties of the liquid: we 
therefore cut off its slowly decaying, large k tail by mul- 
tiplying it in Fourier space with a crossover function Ac- 
given by 



A cr (k) 



1 



(i+ik|/fc c r 



(14) 



with k c — 0.33 bohr, chosen to keep the C functions as 
band-width limited as possible. 

Then we determine C as described above using the 
partial structure factor data for the uniform liquid from 
Figure 1 in [30(. The best, explicit partial structure fac- 
tor data which we have found are measured at standard 
pressure, but at a temperature T = 20 ± 3 °C, which dif- 
fers slightly from the standard temperature (T = 25 °C) 
employed in this work. Comparing experimental radial 
distribution functions at various temperatures [3~0 ]. we 
find that the differences remain smaller than 10% over 
a temperature range of 30 Kelvin. We therefore expect 
an overall error of no more than one or two percent due 
to this difference in temperatures. Next, adjusting the 
smoothing parameter to give the experimental surface 
tension of 71.98 x 10~ 5 N/m 28], measured at standard 
conditions, yields r = 4.2027 bohr (to all the digits 
used in our software). Finally, to provide a computa- 
tionally transferable representation of C, we parametrize 
each component C ai by a sum of Gaussians in Fourier 
space 



C Q7 (k) = ]T AW exp -flW (|k| - C« 



= 1 



(15) 



Table UTT1 summarizes the resulting fitting coefficients (to 
all decimal places used in our software). To reproduce 
the full matrix, note that C is symmetric and that the 
identity of the hydrogens imposes C\\ = Cxi- Figure [T] 



A (i) 
[hartree] 



B, 



07 

bohr 2 l 



[bohr" 1 ] 





1 


-0.0271 1530 


1.0569 0000 


0.0204 8180 




2 


-0.0795 5700 


1.6385 3000 


0.1295 1300 


Coo 


3 


0.0966 4800 


1.8781 5000 


0.0736 7760 


4 


-0.0291 5170 


2.4678 7000 


0.0563 6360 




5 


0.0227 0520 


3.1004 7000 


0.0844 1320 




6 


-0.0109 0780 


3.7162 0000 


0.0452 0230 




1 


-0.0080 1259 


1.1985 9000 


0.0280 2800 




2 


0.0379 8360 


1.7289 7000 


0.1175 3400 


Coi 


3 


-0.0380 7370 


2.1583 9000 


0.2133 3600 


4 


0.0254 5760 


2.7211 2000 


0.1541 5700 




5 


-0.0029 1329 


3.2989 4000 


0.0265 2820 




6 


0.0010 9967 


3.7659 0000 


0.0315 5180 



Cn 


1 

2 


-0.0139 5900 
0.0295 7760 


1.8869 7000 
2.5316 4000 


0.1041 8600 
0.0869 8480 


Cl2 


1 

2 


-0.0139 5900 
0.0295 7760 


1.8869 7000 
2.5316 4000 


0.1041 8600 
0.0869 8480 



TABLE III: Coefficients of Gaussian parametrized C func- 
tions. Note that, by symmetry, C a - t = C la and Cn = Cii. 
(C12 and Cn turned out to be numerically indistinguishable.) 
Coefficients are given to all digits used in our software to en- 
sure numerical reproducibility of our results. The index 
refers to the oxygen site, the indices 1 and 2 refer to the two 
hydrogen sites. 



> 
O 




2 3 4 5 6 
k [1 /angstrom] 



FIG. 1: Comparison of Coo (thick dashed curve), Coi (thick 
dash-dotted curve) and Cn (thick dotted curve) and the cor- 
responding parametrized forms (thin solid curves) given by 
(|15p and Table iHll Note that C12 is not shown, since it would 
be visually indistinguishable from Cn. 



compares our extracted C functions to the resulting nu- 
merical fits, the latter of which are used in all of the 
calculations presented below. 
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Recipe for Evaluating the Free Energy Functional 
for Water 



In sum, to compute the free-energy functional for a 
given set of effective potentials ^ a (r) , one first evaluates 
the intra-molecular contribution via © and HJ). Next, 
one determines the site densities via ([3]), which allows 
one to evaluate the second term in (JXJ) trivially. Finally, 
one evaluates the last term in {T]) according to {3]), with 
K given by (|5|), C parameterized by the coefficients in 
Table [Ell and F ex computed according to ([5]). ([12 ]) . ([15 ]) . 



WATER IN A PARALLEL PLATE CAPACITOR 

To demonstrate the ability of our density-functional 
theory to describe water in an inhomogeneous environ- 
ment and to capture the dielectric response of the liquid, 
we study the behavior of the liquid in a parallel plate ca- 
pacitor. We carry out our study at standard temperature 
and pressure described above, taking the system to be ho- 
mogeneous in the two dimensions parallel to the capacitor 
plates. Along the perpendicular direction, we impose pe- 
riodic boundary conditions, defining a unit cell contain- 
ing two parallel plate capacitors with opposing externally 
applied electric fields. This arrangement makes the elec- 
trostatic potential a periodic function and eliminates un- 
desired electrostatic interactions between the capacitors. 
Each capacitor has a plate separation of 106 A (200 bohr) 
and the total length of the cell is 423 A (800 bohr). 

In our study we consider both purely repulsive, "hard- 
wall" plates as well as "attractive" plates with a region 
of attraction near the plates. For the hard-wall case, the 
plates are essentially infinite potential hard walls act- 
ing on both the oxygen and hydrogen sites of the wa- 
ter molecules. As a practical numerical matter, to re- 
duce aliasing effects in the use of the numerical Fourier 
transform we approximate such walls by purely repul- 
sive Gaussian potentials of sufficiently narrow width and 
large amplitude so that the results below are insensitive 
to the width and amplitude of the Gaussian employed. 
In particular, we use 



>hw(z) = Aexp ( 



(16) 



with A — 150 kgT and a — 0.5 A. For the attractive case, 
we employ a "9-3" (planarly integrated Lennard- Jones) 
potential acting on the oxygens only, specifically, 



t (z) = C 



D 



z + z shift 



D 



Z + Z s utt 



(17) 



withC= 9x0.15/(2\/3) eV, D = ^ A and z shift =2 A. 
This potential has its minimum at a distance of 3 A from 
the wall at a depth 0.15 eV. We choose this depth as 




5 10 15 20 

distance from capacitor plate [angstrom] 

FIG. 2: Hydrogen- (thick dotted curve) and oxygen- (thick 
solid curve) site density in weak field and zero field (thin 
curves) versus distance from hard-wall plate. (Note that the 
hydrogen and oxygen densities are nearly visually indistin- 
guishable.) 



typical of the interactions in such systems. For instance, 
Berkowitz and co-workers [32| employ a corrugated po- 
tential with an average depth of about 0.40 eV, which is 
also in the range of a couple of tenths of an electron Volt. 

To determine the thermodynamic state of the system 
we then expand the density fields on a real space grid with 
8192 sampling points (~ 20 points/A), and minimize the 
grand free energy |lj using standard numerical conjugate 
gradient techniques without any preconditioning. 



Results 

Hard-wall plates 

We begin with a brief discussion of our results for the 
hard wall case as an idealized reference for comparison. 
Figure [21 shows our results for oxygen and hydrogen equi- 
librium density profiles in the absence of an externally 
applied electrostatic displacement (D = 0) and compares 
them to the density profiles of hydrogen and oxygen sites 
in a relatively weak external field (a = (3pD/e ps 0.009). 
The zero-field profiles exhibit an extended gas phase re- 
gion up to a distance of ~ 10 A from the plates. This 
"lingering" gas phase region exists because molecules can, 
at very little free-energy cost, minimize the influence of 
the repulsive plates by leaving the system. (Recall that 
([1]) is written in the grand canonical ensemble, so that 
the system can exchange particles with an external reser- 
voir.) As soon as a relatively weak external electric dis- 
placement is applied, however, it becomes favorable for 
the dipolar molecules to enter the system because they 
can lower their energy by partially aligning with the field 
within the capacitor. For even quite small fields, this ef- 
fect causes the capacitor to fill with liquid, thereby elim- 
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FIG. 3: Hydrogen- (dashed curve) and oxygen- (solid curve) 
site density in strong field versus distance from hard-wall 
plate. 

mating most of the gas phase region. The result of this 
is then an almost rigid shift of the density profiles to- 
ward the capacitor plates, with an appearance and gen- 
eral shape which remains nearly fixed over a large range 
of fields. 

Considering the fine details of the density profiles, we 
note that, contrary to our previous results for diatomic 
liquids, a relatively small (barely noticeable in the figure) 
asymmetry exists at zero external field between the oxy- 
gen and hydrogen site densities so that there is a non- 
vanishing electrostatic charge density in the vicinity of 
the capacitor plates — even in the absence of an exter- 
nally applied electric field. (Note: because the oxygen 
site carries a charge which is double in magnitude com- 
pared to the hydrogen site and because there are two in- 
dependent but identical hydrogen sites, the magnitude of 
the charges carried by the oxygen and hydrogen density 
fields are directly proportional to the curves in the figure, 
with the same constant of proportionality for each.) By 
symmetry, the net charge which the plates induce must 
integrate to zero; however, a finite surface dipole remains 
with a corresponding potential jump which we find in the 
zero-field hard-wall case to be +0.041 V across the solid- 
liquid interface, with the potential being higher in the 
solid. 

Figure [3] shows the equilibrium density profiles in a rel- 
atively strong (a w 0.25) external field directed from left 
to right in the figure. The relative abundance of oxygen 
sites close to the left plate (which carries a positive ex- 
ternal charge to generate the field) and of hydrogen sites 
next to the right plate (which carries a negative exter- 
nal charge) indicates that the water molecules exhibit a 
strong tendency to align their dipoles with the applied 
field. Interestingly, this leads to quite different density 
profiles on the opposite sides of the capacitor. On the 
left (positive) plate, the first peak in the oxygen density 
is followed by a peak in the hydrogen density of similar 
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FIG. 4: Hydrogen- (thick dashed curve) and oxygen- (thick 
solid curve) site density in strong field versus distance from 
attractive wall. Reference result for the zero field oxygen 
density at repulsive wall (thin solid curve). 

shape, width and height. This we interpret to be a di- 
rect consequence of alignment of the molecules with the 
field and the rigid intra-molecular bonding between hy- 
drogen and oxygen atoms, an effect quite similar to what 
was observed in our previous study of liquid hydrogen 
chloride [23|]. The distance between these peaks, 0.65 A, 
is within about 10% of 0.58 A, what one expects from 
the basic geometry of the water molecule assuming that 
the molecules fully align with the field. 

On the right (negative) plate, we again find a well- 
defined peak in the oxygen density, though the peak is 
somewhat stronger than on the left. We find also that the 
hydrogen density, rather than being also sharply peaked, 
is now spread significantly around the oxygen peak. This 
case is different, because the hydrogen sites, now pointing 
away from the liquid region, cannot form hydrogen bonds 
if the molecules align fully with the field. 



Attractive plates 

Figure [4] compares the zero-field results for the hard- 
wall plates with those of the attractive plates. There 
is no longer a "lingering" gas phase region because the 
minimum in the interaction potential (acting only on the 
oxygen sites in this case) with the plates sets the loca- 
tion of the main peak in the oxygen density. This mini- 
mum, as described above, occurs at a distance of about 
3 A from the attractive walls and corresponds well with 
the peaks in the oxygen and hydrogen density profiles. 
The height of the oxygen peak is larger than the den- 
sity in the bulk by a factor of ~ 4.7, which is similar 
to the results of molecular dynamics calculations of wa- 
ter adjacent to platinum surfaces [32| and experimental 
findings [33j ]. Finally, there is also a net dipole layer in 
our attractive case leading now to a significantly stronger 
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FIG. 5: Hydrogen- (dashed curve) and oxygen- (solid curve) 
site density in strong field versus distance from attractive wall. 



potential jump at the liquid-solid interface of +0.220 V, 
with the potential in the solid again being higher. This 
value is typical of realistic systems: from molecular dy- 
namics data of water adjacent to a platinum surface 34fl . 
we extract a potential jump of ~ +0.5 V with the poten- 
tial being higher in the solid. 

When a strong field (again, a w 0.25) is applied to the 
capacitor with attractive walls (Figure [5]), we find a rear- 
rangement of the liquid structure which exhibits various 
quantitative differences compared to the hard-wall case: 
the oxygen peaks are now significantly more pronounced 
(reaching values near 4 and 5 times the equilibrium liquid 
density, as opposed to the previous peak values near 1.5 
and 2.2) as the oxygen atoms settle into the attractive 
potential wells near the walls. Also, the peak near the 
positive plate is also now notably broader (with a width 
of 1.6 A at the value of the equilibrium, as opposed to 
the previous value of 1.0 A). Finally, the hydrogen peak is 
significantly less pronounced relative to the oxygen peak 
and the peak position is now less displaced (0.4 A) from 
the oxygen peak position. Despite these quantitative dif- 
ferences, however, the qualitative features at strong fields 
are quite similar in the attractive and the hard- wall cases. 
Namely, on the left (positive) plate we find a well-defined 
peak in the oxygen density followed by a corresponding 
peak in the hydrogen density, and on the right (negative) 
plate there is a well-defined oxygen peak which is sur- 
rounded by a more diffuse peak in the hydrogen density. 
Also, as with the hard wall case, we find that the hydro- 
gen peak is off-center from the oxygen peak and displaced 
toward the negative plate to which it is attracted. 

This suggests that, at these large field strengths, gen- 
eral features are independent of the exact form of the 
liquid-solid interaction. In fact, we note that this sort 
of asymmetry between positive and negative plates is 
also evident in molecular dynamics calculations which, 
at fields near a = 0.25, also tend to show that, near the 
positive plate, the oxygen density exhibits peaks which 
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FIG. 6: Dielectric function versus total field: density- 
functional results (crosses), non-linear electrostatics (solid 
curve), analytical result of [3(j (dashed line), molecu- 
lar dynamics results of [U (circles), [13] (squares) and 
[3^ | (triangles) . 



arc followed by peaks of similar appearance in the hydro- 
gen density, whereas, near the negative plate, there is a 
well-defined peak in the oxygen density but that the hy- 
drogen density exhibits either a broader peak or multiple 
sub-peaks in the vicinity of the main oxygen peak[34T[35|. 



Nonlinear dielectric response 

The equilibrium site densities also determine the di- 
electric function of water e(D) = D/E(D) with E being 
the total electric field given by the sum of the electric 
displacement and the induced electric field. Although 
E has a complicated form close to the walls, it rapidly 
approaches the constant value E = D — AttP as the dis- 
tance from the walls increases, where P is the induced 
polarization. In the present one-dimensional geometry, 
P is proportional to the induced surface charge which 
appears in the equilibrium density profiles. 

Figure [5] shows our calculated results for e plotted 
as a function of E in units of MV/m (equivalent to 
1 mV/nm or 0.1 mV/A) to simplify comparison with 
the results of previous studies. The figure shows that 
the results from our density-functional theory are well- 
described by self-consistcntly screened nonlinear elec- 
trostatics (solid curve), which gives the polarization P 
as the self-consistent solution to the equation P = 
PW(I) - a e 47rP), where P^ m \E) is the response of 
a gas of noninteracting dipoles in a local field E and 
a e = e/(e - 1) - e ( ni )/( e ( ni ) - 1) ensures that the cor- 
rect linear response is recovered when D is weak. The 
solid curve reproduces our data well and so represents a 
convenient analytic guide to the eye for our results. 

Figure [5] compares our results to Booth's analytical re- 
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suit for the nonlinear dielectric response of water [36j |. 
This result (clashed curve in the figure) was obtained by 
extending Onsager's and Kirkwood's theories of polar di- 
electrics to high field strengths [36| and also reproduces 
quite well explicit molecular dynamics results for high 
fields. It fails, however, at small fields, eventually giv- 
ing an incorrect linear dielectric constant of about ~ 65. 
We find very good agreement with this established re- 
sult for fields larger than 300 MV/m. We also compare 
our results to various molecular dynamics calculations 
3j|, [3j], [38| and find very good agreement, except at zero 
field, where the results of the molecular dynamics calcu- 
lations often show a wide range of results, some differing 
significantly from the experimental dielectric constant of 
~ 80 depending on the model potentials and numerical 
methods employed (39|. The zero-field dielectric constant, 
we, of course, reproduce exactly by construction. Our re- 
sults, therefore, appear to be reliable over a broad range 
of applied fields. 



SUMMARY AND CONCLUSIONS 

This work demonstrates how to apply a new, gen- 
eral Kohn-Sham classical density-functional approach for 
molecular liquids to what is arguably the liquid of great- 
est scientific interest, water. The resulting theory has the 
quite tractable computational cost of the problem of a 
nonintcracting gas of molecules in a self-consistent exter- 
nal potential and gives an exact account of the entropy 
associated with the alignment of the molecular dipole 
moment with an external field. 

After constructing the functional and giving a com- 
plete numerical parameterization of all quantities so that 
our results are numerically reproducible, we demonstrate 
the tractability of the approach by applying it to the 
behavior of water in parallel plate capacitors, with ei- 
ther hard-wall or attractive plates, over a range of ap- 
plied fields through which the dielectric constant e of wa- 
ter varies by over a factor of two. We provide results 
for the distribution of oxygen and hydrogen atoms near 
the capacitor plates and find an asymmetry between the 
response at negatively and positively charged plates in 
general qualitative agreement with molecular dynamics 
results. From the resulting charge distributions, we ex- 
tract predictions of the dielectric constant of water over 
a wide range of field strengths (0 to 800 MV/m) and find 
results in very good agreement with both previous theo- 
retical and molecular dynamics calculations. The afore- 
mentioned exact treatment of the entropy cost of aligning 
the molecules with the field is critical in producing this 
accurate treatment of dielectric saturation effects. 

With a successful framework in place, there is now 
a firm basis for future improvements. Given that we 
now have a good treatment of the internal structure of 
the water molecule, the primary area remaining for im- 



provement is our choice to employ for this first work a 
somewhat simplistic weighted-density functional form in 
((6|). The direct route to improve upon this is to leverage 
much of the existent work in simple fluids, both in tra- 
ditional weighted density- functional theories [2^] and in 
fundamental measure theory [i(J. Another clear area for 
improvement would be the inclusion of additional, ionic 
species into the functional to allow for screening by elec- 
trolytes, which should be relatively straightforward. 

In sum, we present a new framework for a computa- 
tionally tractable continuum description of water which 
captures, within a unified approach with a firm theo- 
retical foundation, molecular-scale correlations, entropic 
effects, microscopic dielectric screening, surface tension, 
bulk thermodynamic properties, and dielectric satura- 
tion. This description thus includes a priori all of the 
various effects which are generally added a posteriori to 
continuum models of water. Although less detailed than 
current explicit molecular treatments of the solvent, freed 
from the need for computing thermodynamic averages, 
a sufficiently accurate continuum approach such as we 
present has the potential to open up much larger and 
complex systems to theoretical study. As such, the new 
framework is an ideal starting point for a variety of new 
investigations into a wide class of subjects, which include 
protein folding, molecular motors, membrane physics, 
drug design, the electrochemical interface, liquid-phase 
catalysis and fuel cells. 
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